Proyecto final
Alejandra Lelo de Larrea Ibarra, C.U. 000124433
Luis Fernando Cantú Díaz de León, C.U. 000125755
15 de diciembre de 2017
Problema y Objetivo
Los sistemas para compartir bicicletas (como EcoBici en México) son una nueva generación del tradicional alquiler de bicicletas. En éstos, todo el proceso de membresía, alquiler y devolución se ha automatizado. Hoy en día, existe un gran interés en estos sistemas debido a su importante papel en la reducción del tráfico vehicular, el cuidado del medio ambiente y la reducción en los problemas de salud. Actualmente, hay alrededor de 500 programas para compartir bicicletas en todo el mundo que se componen por más de 500 mil bicicletas.
A través de estos sistemas, el usuario puede alquilar una bicicleta fácilmente desde una estación particular y regresarla en una estación que no tiene porqué coincidir con la inicial. Opuesto a otros servicios de transporte como el autobús o el metro, la duración del viaje, la posición de salida y la de llegada se registran explícitamente en estos sistemas. Esta característica convierte el sistema de intercambio de bicicletas en una red de sensores virtuales que se puede usar para detectar la movilidad en la ciudad y otras características. Por ejemplo, el proceso de alquiler de bicicletas compartidas está altamente correlacionado con las condiciones ambientales. Las condiciones climáticas, la precipitación, el día de la semana, la estación del año, la hora del día, etc. pueden afectar las conductas de alquiler.
De esta manera, se busca predecir el número de bicicletas rentadas por hora según las condiciones climáticas y ambientales y comparar los resultados obtenidos con distintos métodos vistos en clase. Para ello, se utiliza el registro histórico (años 2011 y 2012) del sistema Capital Bikeshare, Washington D.C., EE. UU.[^1] La base de datos contiene 17,379 observaciones para 17 variables; entre ellas se encuentran: la fecha y hora, estación del año, si el día es laborable o no, tipo de clima, condiciones ambientales (temperatura, humedad, etc) y número de bicicletas rentadas (variable de interés). Los métodos a utilizar para la estimación y comparación son: regresión lineal, regresión lineal con regularización, redes neuronales, bosques aleatorios y gradient boosting machine. [^1]: Los datos fueron recuperados de UCI Machine Learning Repository.
Datos y análisis descriptivo
# Se leen los datos
datos <- read_csv("hour.csv",
col_types =
cols(
instant = col_integer(),
dteday = col_date(format = ""),
season = col_character(),
yr = col_integer(),
mnth = col_character(),
hr = col_character(),
holiday = col_integer(),
weekday = col_character(),
workingday = col_integer(),
weathersit = col_character(),
temp = col_double(),
atemp = col_double(),
hum = col_double(),
windspeed = col_double(),
casual = col_integer(),
registered = col_integer(),
cnt = col_integer()))
datos_dummies <- dummy.data.frame(datos%>%data.frame, sep = "_") %>% as.tibble() %>%
select(-instant, -casual, -registered)
dim(datos)## [1] 17379 17Variables
Las variables contenidas en la base de datos son:
- instant: índice de observaciones
- dteday : fecha
- season : Estación del año (1=primavera, 2=verano, 3=otoño, 4=invierno)
- yr : año (0=2011, 1=2012)
- mnth : mes ( 1 to 12)
- hr : hora (0 to 23)
- holiday : día festivo (1) o no festivo (0)
- weekday : día de la semana (0 a 6 con 0=Domingo)
- workingday : día laboral (1= no es fin de semana ni festivo, 0 e.o.c.)
weathersit : clima
+ 1 = Claro, Pocas nubes, Parcialmente nublado, Parcialmente nublado + 2 = Niebla + Nublado, Niebla + Nubes fragmentadas, Niebla + Pocas nubes, Niebla + 3 = Nieve ligera, lluvia ligera + tormenta eléctrica + nubes dispersas, lluvia ligera + nubes dispersasred clouds + 4 = luvia pesada + paletas de hielo + tormenta + niebla, nieve + nieblaist, Snow + Fog- temp : Temperatura normalizada en grados Celcius. Divididos entre 41 (máximo)
- atemp: Sensación térmica normalizada en grados Celcius. Divididos entre 50 (máximo)
- hum: Humedad normalizada. Divididos entre 100 (máximo)
- windspeed: velocidad del viendo normalizada. Divididos entre 67 (máximo)
- casual: número de usuarios casuales
- registered: número de usuarios registrados
cnt: número total de bicicletas rentadas (incluye usuarios casuales y registrados)
Manipulación de datos
En primer lugar, eliminamos las variables casual y registered pues son variables que no se tendrán en la verdadera tarea de predicción. Transformarmos las variables categóricas (season, mnth, hr, weekday y weathersit ) en dummies para facilitar su manejo en algunos modelos.
# Conservamos la base con datos sin one hot encoding
datos_sin <- datos %>%
mutate(datetime = as.POSIXct((dteday) + hours(hr))) %>%
select(-casual, -registered, -instant)
# Eliminamos las variables casual y registered, también creamos la variable datetime, para poder hacer una mejor visualización de los datos.
datos<-datos%>%
select(-c(casual,registered))
# Almacenamos la fecha
fecha<-paste(datos$dteday,datos$hr,sep="/")
# Se asignan como factores las variables casual y registred pues
cols<-match(c('season', 'mnth', 'hr', 'weekday','weathersit'),names(datos))
datos[cols]<-lapply(datos[,cols],as.factor)
# Se crea una dummy para cada categoría (menos una)
auxdatos<-model.matrix(~0+season+mnth+hr+weekday+weathersit,datos)
# Se eliminan las varaibles categóricas y se agregan las dummies
datos<-datos%>%
select(-c(season, mnth, hr, weekday,weathersit))%>%
mutate(dteday = ymd(dteday)) %>%
cbind(auxdatos)
head(datos) ## instant dteday yr holiday workingday temp atemp hum windspeed cnt
## 1 1 2011-01-01 0 0 0 0.24 0.2879 0.81 0.0000 16
## 2 2 2011-01-01 0 0 0 0.22 0.2727 0.80 0.0000 40
## 3 3 2011-01-01 0 0 0 0.22 0.2727 0.80 0.0000 32
## 4 4 2011-01-01 0 0 0 0.24 0.2879 0.75 0.0000 13
## 5 5 2011-01-01 0 0 0 0.24 0.2879 0.75 0.0000 1
## 6 6 2011-01-01 0 0 0 0.24 0.2576 0.75 0.0896 1
## season1 season2 season3 season4 mnth10 mnth11 mnth12 mnth2 mnth3 mnth4
## 1 1 0 0 0 0 0 0 0 0 0
## 2 1 0 0 0 0 0 0 0 0 0
## 3 1 0 0 0 0 0 0 0 0 0
## 4 1 0 0 0 0 0 0 0 0 0
## 5 1 0 0 0 0 0 0 0 0 0
## 6 1 0 0 0 0 0 0 0 0 0
## mnth5 mnth6 mnth7 mnth8 mnth9 hr1 hr10 hr11 hr12 hr13 hr14 hr15 hr16
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 1 0 0 0 0 0 0 0
## 3 0 0 0 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0 0 0 0 0
## hr17 hr18 hr19 hr2 hr20 hr21 hr22 hr23 hr3 hr4 hr5 hr6 hr7 hr8 hr9
## 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
## 3 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0
## 4 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
## 5 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0
## 6 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0
## weekday1 weekday2 weekday3 weekday4 weekday5 weekday6 weathersit2
## 1 0 0 0 0 0 1 0
## 2 0 0 0 0 0 1 0
## 3 0 0 0 0 0 1 0
## 4 0 0 0 0 0 1 0
## 5 0 0 0 0 0 1 0
## 6 0 0 0 0 0 1 1
## weathersit3 weathersit4
## 1 0 0
## 2 0 0
## 3 0 0
## 4 0 0
## 5 0 0
## 6 0 0Separación de muestras
La separación de muestras se realizará sobre la línea de tiempo: 1. Muestra de entrenamiento: datos por hora del 1° de enero del 2011 al 30 de septiembre del 2012. 2. Muestra de validación: datos por hora del 1° de octubre del 2012 al 30 de noviembre del 2012. 3. Muestra de prueba: datos por hora del 1 de diciembre del 2012 al 31 de diciembre del 2012.
test <- datos %>% filter(dteday >= "2012-12-01") %>% select(-instant, -dteday)
valid <- datos %>% filter(dteday >= "2012-10-01" & dteday < "2012-11-30") %>%
select(-instant, -dteday)
train <- datos %>% filter(dteday < "2012-09-30") %>% select(-instant, -dteday)
test_dummies <- datos_dummies %>% filter(dteday >= "2012-12-01") %>% select(-dteday)
valid_dummies <- datos_dummies %>% filter(dteday >= "2012-10-01" & dteday < "2012-11-30") %>%
select(-dteday)
train_dummies <- datos_dummies %>% filter(dteday < "2012-09-30") %>%
select(-dteday)
# Datos sin one hot encoding
test_sin <- datos_sin %>% filter(dteday >= "2012-12-01") %>% select(-datetime, -dteday)
valid_sin <- datos_sin %>% filter(dteday >= "2012-10-01" & dteday < "2012-11-30") %>%
select(-datetime, -dteday)
train_sin <- datos_sin %>% filter(dteday < "2012-09-30") %>% select(-datetime, -dteday)
# Fechas almacenadas por aparte
fechas_test <- datos_sin %>% filter(dteday >= "2012-12-01") %>% select(datetime, dteday)
fechas_valid <- datos_sin %>% filter(dteday >= "2012-10-01" & dteday < "2012-11-30") %>%
select(datetime, dteday)
fechas_train <- datos_sin %>% filter(dteday < "2012-09-30") %>% select(datetime, dteday)Una vez separadas las muestras, se estandarizan las variables numéricas (menos la respuesta) en los datos de entrenamiento y, con esas medias y varianzas, se estandarizan los datos de validación y prueba.
Visualización de datos
Primero podemos ver el número de usuarios de las bicicletas por cada hora, tomando en cuenta la temperatura y la estación del año en que se utiliza. Se puede apreciar que, como podría esperarse, hay un mayor número de usuarios entre las 7:00 y 8:00 AM, que es a la hora en que las personas se trasladan a su trabajo, y también entre 5:00 y 6:00 PM, cuando van de regreso a sus hogares. Además se puede ver que se utilizan menos las bicicletas durante las estaciones más frías, invierno y primavera.
graf <- ggplot(data = cbind(train_sin, fechas_train) )
graf +
geom_jitter(aes(x = as.numeric(hr), y = cnt, color = temp)) +
ggtitle("Número de usuarios de acuerdo a la hora") +
xlab("Hora") +
ylab("Número de usuarios") +
scale_x_continuous(breaks=c(0:23), labels=c(0:23)) +
facet_grid(season ~ .,
labeller = as_labeller(c("1" = "Primavera","2" = "Verano",
"3" = "Otoño","4" = "Invierno"))) +
theme_tufte(base_size = 14, base_family = "sans")
La estación del año podría ser de las variables importantes a la hora de que un usuario decida utilizar una bicicleta o no, ya que esperaríamos que los usuarios utilicen el servicio en las temporadas más cálidas. Esto lo podemos comprobar en la siguiente gráfica, donde se puede ver que hay un promedio de usuarios un poco más bajo durante la primavera y el invierno. Además de eso, también se aprecia un promedio menor los sábados y domingos.
ggplot() +
geom_boxplot(data = rbind(train_sin, valid_sin), aes(x = season, y = cnt, fill =as.factor(weekday))) +
ggtitle("Número de usuarios de acuerdo a la estación del año") +
xlab("Estación del año") +
ylab("Número de usuarios") +
scale_x_discrete(labels = c("Primavera", "Verano", "Otoño", "Invierno")) +
scale_fill_discrete(name = "Día",
labels = c("Domingo", "Lunes", "Martes",
"Miércoles", "Jueves", "Viernes", "Sábado")) +
theme_tufte(base_size = 14, base_family = "sans")
También las condiciones de lluvia y humedad pueden afectar el número de personas que deciden transportarse en bicicleta. En la siguiente gráfica se puede apreciar que esta relación sí existe para los días laborales: a mayor humedad, menor número de usuarios.
graf +
geom_jitter(aes(x = as.numeric(hr), y = cnt, color = hum)) +
ggtitle("Número de usuarios y humedad") +
xlab("Hora") +
ylab("Número de usuarios") +
scale_x_continuous(breaks=c(0:23), labels=c(0:23)) +
facet_grid(workingday~.,
labeller = as_labeller(c("1" = "Laboral","0" = "No laboral"))) +
theme_tufte(base_size = 14, base_family = "sans")
Ajuste y Validación de modelos
En esta sección realizamos el ajuste del modelo con distintos métodos vistos en clase. Dentro de cada método, la selección de hiperparámetros se realiza con base en el menor error de validación.
Regresión Lineal
Para realizar la regresión lineal, usamos los datos sin one hot encoding.
modelo_lineal <- lm(cnt ~ . -workingday, data = train_sin)
summary(modelo_lineal)##
## Call:
## lm(formula = cnt ~ . - workingday, data = train_sin)
##
## Residuals:
## Min 1Q Median 3Q Max
## -354.14 -59.34 -6.58 49.72 437.50
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -72.296 6.794 -10.641 < 2e-16 ***
## season2 29.753 4.860 6.122 9.48e-10 ***
## season3 16.615 6.061 2.741 0.006128 **
## season4 42.214 6.029 7.002 2.63e-12 ***
## yr 89.030 1.767 50.393 < 2e-16 ***
## mnth10 42.722 7.958 5.369 8.05e-08 ***
## mnth11 28.400 7.808 3.637 0.000276 ***
## mnth12 29.368 6.107 4.809 1.53e-06 ***
## mnth2 5.477 3.818 1.434 0.151451
## mnth3 22.142 4.356 5.083 3.77e-07 ***
## mnth4 21.640 6.566 3.296 0.000983 ***
## mnth5 39.269 7.095 5.535 3.17e-08 ***
## mnth6 29.884 7.430 4.022 5.79e-05 ***
## mnth7 16.840 8.508 1.979 0.047807 *
## mnth8 36.507 8.284 4.407 1.05e-05 ***
## mnth9 61.087 7.579 8.060 8.20e-16 ***
## hr1 -17.507 5.559 -3.149 0.001639 **
## hr10 105.334 5.585 18.859 < 2e-16 ***
## hr11 129.542 5.626 23.026 < 2e-16 ***
## hr12 166.760 5.673 29.397 < 2e-16 ***
## hr13 162.868 5.716 28.492 < 2e-16 ***
## hr14 147.541 5.747 25.672 < 2e-16 ***
## hr15 155.252 5.761 26.949 < 2e-16 ***
## hr16 215.865 5.749 37.545 < 2e-16 ***
## hr17 370.634 5.720 64.801 < 2e-16 ***
## hr18 343.219 5.684 60.385 < 2e-16 ***
## hr19 238.309 5.627 42.349 < 2e-16 ***
## hr2 -25.709 5.581 -4.607 4.13e-06 ***
## hr20 158.771 5.594 28.383 < 2e-16 ***
## hr21 109.864 5.570 19.724 < 2e-16 ***
## hr22 72.506 5.558 13.045 < 2e-16 ***
## hr23 32.542 5.554 5.860 4.74e-09 ***
## hr3 -36.799 5.620 -6.548 6.01e-11 ***
## hr4 -39.986 5.630 -7.102 1.28e-12 ***
## hr5 -23.701 5.590 -4.240 2.25e-05 ***
## hr6 33.919 5.573 6.087 1.18e-09 ***
## hr7 164.757 5.561 29.628 < 2e-16 ***
## hr8 298.289 5.554 53.703 < 2e-16 ***
## hr9 156.913 5.561 28.217 < 2e-16 ***
## holiday -21.005 5.292 -3.970 7.24e-05 ***
## weekday1 4.985 3.103 1.606 0.108210
## weekday2 7.819 3.021 2.588 0.009651 **
## weekday3 9.348 3.024 3.091 0.001998 **
## weekday4 9.698 3.022 3.209 0.001335 **
## weekday5 13.034 3.015 4.323 1.55e-05 ***
## weekday6 12.573 2.999 4.193 2.77e-05 ***
## weathersit2 -9.653 2.031 -4.753 2.02e-06 ***
## weathersit3 -61.611 3.336 -18.469 < 2e-16 ***
## weathersit4 -62.670 57.265 -1.094 0.273804
## temp 124.989 29.747 4.202 2.66e-05 ***
## atemp 87.224 30.725 2.839 0.004534 **
## hum -82.439 5.694 -14.477 < 2e-16 ***
## windspeed -31.182 7.290 -4.277 1.90e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 98.9 on 15134 degrees of freedom
## Multiple R-squared: 0.6915, Adjusted R-squared: 0.6904
## F-statistic: 652.2 on 52 and 15134 DF, p-value: < 2.2e-16Se calcula el error cuadrático medio del modelo lineal.
preds_lineal <- predict(modelo_lineal, newdata = valid_sin)
resid_lineal <- valid_sin$cnt - preds_lineal
mse_valid_lineal <- mean(resid_lineal**2)
mse_valid_lineal## [1] 16486.67Se elaboran gráficas para observar el desempeño de este modelo. En la primera gráfica se puede apreciar que el modelo generalmente identifica la tendencia temporal, pero tiende a subestimar el número real de usuarios que utilizan las bicicletas. En la segunda gráfica podemos observar que los residuos del modelo se alejan mucho del 0, confirmando la subestimación. En general, se puede observar que no se logra un buen ajuste de los datos de validación.
ggplot() +
geom_line(aes(x = fechas_valid$datetime, y = valid_sin$cnt), color = "blue", alpha = 0.5, size = 0.5) +
geom_line(aes(x = fechas_valid$datetime, y = preds_lineal), color = "red", alpha = 0.5, size = 0.5) +
ggtitle("Predicciones del modelo lineal y valores reales a través del tiempo") +
ylab("Número de usuarios") +
xlab("Fecha") +
theme_tufte(base_size = 14, base_family = "sans")
ggplot() +
geom_point(aes(x = fechas_valid$datetime, y = resid_lineal)) +
geom_abline(intercept = 0, slope = 0, color = "red", size = 1.2) +
ggtitle("Residuales del modelo lineal") +
ylab("Residuales") +
xlab("Fecha") +
theme_tufte(base_size = 14, base_family = "sans")
ggplot() +
geom_jitter(aes(x = valid$cnt, y = preds_lineal)) +
geom_abline(intercept = 0, slope = 1, color = "red", size = 1.2) +
labs(x = "Usuarios", y = "Usuarios predecidos") +
ggtitle("Predicciones del modelo lineal vs. Valores reales") +
theme_tufte(base_size = 14, base_family = "sans")
Regresión Lineal con regularización
Creamos una función que ajuste un nuevo modelo de acuerdo a cada lambda que se encuentra dentro de una lista que se utiliza como input.
modelos_lasso <- function(lambdas, alphas){
dfs <- data.frame()
x_train <- train_dummies %>% select(-cnt) %>% data.matrix()
y_train <- train_dummies %>% select(cnt) %>% data.matrix()
x_valid <- valid_dummies %>% select(-cnt) %>% data.matrix()
y_valid <- valid_dummies %>% select(cnt) %>% data.matrix()
x_test <- test_dummies %>% select(-cnt) %>% data.matrix()
y_test <- test_dummies %>% select(cnt) %>% data.matrix()
for(j in alphas){
for(i in lambdas){
set.seed(576327)
mod_glm <- glmnet(x = x_train,
y = y_train,
lambda = i,
family = "gaussian",
intercept = T,
alpha = j)
# Calculamos los residuos
preds_train <- predict(mod_glm, newx = x_train)
resid_train <- y_train - preds_train
preds_valid <- predict(mod_glm, newx = x_valid)
resid_valid <- y_valid - preds_valid
# Calculamos el error cuadrático medio
mse_train <- mean((resid_train)**2) %>% sqrt()
mse_valid <- mean((resid_valid)**2) %>% sqrt()
# Se crea una dataframe con info de cada modelo
df <- data.frame(lambdas = i,
alphas = j,
MSE_Train = mse_train,
MSE_Valid = mse_valid
)
# Agregamos el dataframe anterior a dfs
dfs <- rbind(dfs, df)
}
}
# Tomamos la lambda que minimiza el error cuadrático medio de valid
lambda_min <- dfs$lambdas[which.min(dfs$MSE_Valid)]
alpha_min <- dfs$alphas[which.min(dfs$MSE_Valid)]
# Corremos el modelo con los hiperparámetros anteriores
glm_mejor <- glmnet(x = x_train,
y = y_train,
lambda = lambda_min,
family = "gaussian",
intercept = T,
standardize = F,
alpha = alpha_min)
# Generamos los residuos de validación del mejor modelo
preds <- predict(glm_mejor, newx = x_valid)
resids <- y_valid - preds
return(list(modelos = dfs,
mejor = glm_mejor,
x_train = x_train,
x_valid = x_valid,
x_test = x_test,
y_train = y_train,
y_valid = y_valid,
y_test = y_test,
preds = preds,
resids = resids))
}Corremos 222 modelos, tanto LASSO como ridge, y los ordenamos por error cuadrático medio de validación.
modelos_lass <- modelos_lasso(lambdas = exp(seq(-10,10,.05)), alphas = c(0,1))
modelos_lass[[1]] %>% arrange(MSE_Valid) %>% head(10) %>% knitr::kable()| lambdas | alphas | MSE_Train | MSE_Valid |
|---|---|---|---|
| 4.54e-05 | 0 | 98.73236 | 128.3908 |
| 4.77e-05 | 0 | 98.73236 | 128.3908 |
| 5.02e-05 | 0 | 98.73236 | 128.3908 |
| 5.27e-05 | 0 | 98.73236 | 128.3908 |
| 5.55e-05 | 0 | 98.73236 | 128.3908 |
| 5.83e-05 | 0 | 98.73236 | 128.3908 |
| 6.13e-05 | 0 | 98.73236 | 128.3908 |
| 6.44e-05 | 0 | 98.73236 | 128.3908 |
| 6.77e-05 | 0 | 98.73236 | 128.3908 |
| 7.12e-05 | 0 | 98.73236 | 128.3908 |
Graficamos las lambdas para cualesquier alpha. La línea roja vertical indica la lambda con la que se minimiza el error cuadrático medio de validación.
ggplot() +
facet_grid(alphas~ .) +
geom_point(data = modelos_lass[[1]], aes(x = log(lambdas), y = (MSE_Valid), color =alphas), size=0.7) +
geom_vline(xintercept = modelos_lass[[1]]$lambda[which.min(modelos_lass[[1]]$MSE_Valid)] %>% log(),
color = "red") +
ylab("MSE Validación") +
theme_tufte(base_size = 14, base_family = "sans")
La lambda que minimiza el error cuadrático medio de validación es muy pequeña, por lo que los resultados de este modelo se parecen mucho a los de la regresión lineal de la sección anterior.
ggplot() +
geom_line(aes(x = fechas_valid$datetime, y = modelos_lass$y_valid), color = "blue") +
geom_line(aes(x = fechas_valid$datetime, y = modelos_lass$preds), color = "red") +
ylab("Número de usuarios") +
xlab("Fecha") +
theme_tufte(base_size = 14, base_family = "sans")
ggplot() +
geom_point(aes(x = fechas_valid$datetime, y = modelos_lass$resids)) +
geom_abline(intercept = 0, slope = 0, color = "red", size = 1.2) +
ggtitle("Residuales del modelo LASSO") +
ylab("Residuales") +
xlab("Fecha") +
theme_tufte(base_size = 14, base_family = "sans")
ggplot() +
geom_point(aes(x = valid$cnt, y = modelos_lass$preds)) +
geom_abline(intercept = 0, slope = 1, color = "red", size = 1.2) +
labs(x = "Usuarios", y = "Usuarios predecidos") +
ggtitle("Predicciones de LASSO vs. Valores reales") +
theme_tufte(base_size = 14, base_family = "sans")
Redes Neuronales
Ahora ajustamos un modelo de redes neuronales.
modelos_redes <- function(n.units, dropout_rate = 0){
x_train <- modelos_lass$x_train
y_train <- modelos_lass$y_train
x_valid <- modelos_lass$x_valid
y_valid <- modelos_lass$y_valid
models <- list()
history <- list()
for(j in dropout_rate){
for(i in n.units){
set.seed(576327)
model <- keras_model_sequential()
model %>%
layer_dense(units = i,
activation = 'relu',
kernel_initializer=initializer_random_uniform(minval=-0.5,maxval=0.5),
input_shape = c(ncol(x_train))) %>%
layer_dropout(rate = j) %>%
layer_dense(units = floor(i/2),
activation = 'relu',
kernel_initializer=initializer_random_uniform(minval=-0.5,maxval=0.5)) %>%
layer_dropout(rate = j) %>%
layer_dense(units = 1)
model %>% compile(
loss = 'mse',
optimizer = "adam",
metrics = c('mse','mae')
)
iteraciones <- model %>% fit(
x_train,
y_train,
epochs = 1000,
batch_size= nrow(x_train),
verbose = 0,
validation_data = list(x_valid,y_valid)
)
models <- c(models, model)
history <- c(history, iteraciones)
}
}
return(list(models = models,
histories = history))
}Calculamos los residuos de la red neuronal.
modelos_neuronales <- modelos_redes(c(124, 256, 512), c(0, 0.3, 0.5))
preds_neuronales <- modelos_neuronales$models[[9]] %>% predict_on_batch(modelos_lass$x_valid)
resid_neuronales <- modelos_lass$y_valid - preds_neuronalesGraficamos los resultados anteriores. La red neuronal parece dar mucho mejores resultados que los modelos intentados previamente, aunque en la segunda gráfica se puede apreciar que este modelo no parece generalizar correctamente algunas observaciones entre el 15 y el 30 de noviembre.
ggplot() +
geom_line(aes(x = fechas_valid$datetime, y = modelos_lass$y_valid), color = "blue", alpha = 0.5,
size = 0.5) +
geom_line(aes(x = fechas_valid$datetime, y = preds_neuronales), color = "red", alpha = 0.5, size = 0.5) +
ggtitle("Predicciones de red neuronal y valores reales a través del tiempo") +
ylab("Número de usuarios") +
xlab("Fecha") +
theme_tufte(base_size = 14, base_family = "sans")
ggplot() +
geom_point(aes(x = fechas_valid$datetime, y = resid_neuronales)) +
geom_abline(intercept = 0, slope = 0, color = "red", size = 1.2) +
ggtitle("Residuales de red neuronal") +
ylab("Residuales") +
xlab("Fecha") +
theme_tufte(base_size = 14, base_family = "sans")
ggplot() +
geom_point(aes(x = valid$cnt, y = preds_neuronales)) +
geom_abline(intercept = 0, slope = 1, color = "red", size = 1.2) +
labs(x = "Usuarios", y = "Usuarios predecidos") +
ggtitle("Predicciones de red neuronal vs. Valores reales") +
theme_tufte(base_size = 14, base_family = "sans")
Bosque Aleatorio
Construimos una función para correr modelos con diferentes parámetros:
bosques <- function(train, valid, ntrees){
mtry <- c(1:12)
forest <- list()
dfs <- data.frame()
for(i in mtry){
set.seed(576327)
bosque <- randomForest(cnt ~ .,
data = train_sin,
ntree = ntrees,
mtry = i,
importance=TRUE)
# Encontramos las predicciones que arroja el modelo
# probs_train <- predict(bosque, newdata = train)
probs_valid <- predict(bosque, newdata = valid_sin)
# Calculamos el error cuadrático medio
mse_valid <- ((valid_sin$cnt - probs_valid)^2) %>% mean()
df <- data.frame(mtry = i,
ntrees = ntrees,
MSE_Valid = mse_valid
)
dfs <- rbind(dfs, df)
forest <- c(forest, list(bosque))
}
# return(forest)
return(dfs)
}number_trees <- function(ntrees = c(100, 500)){
dfs <- data.frame()
for(i in ntrees){
forest <- bosques(train, valid, i)
dfs <- rbind(dfs, forest)
}
return(dfs)
}Corremos los modelos:
modelos_forests <- number_trees()Los ordenamos de acuerdo al error cuadrático medio de validación.
modelos_forests %>% arrange(MSE_Valid) %>% head(10) %>% knitr::kable()| mtry | ntrees | MSE_Valid |
|---|---|---|
| 10 | 500 | 5116.339 |
| 9 | 500 | 5174.625 |
| 9 | 100 | 5189.864 |
| 11 | 500 | 5193.212 |
| 8 | 100 | 5199.832 |
| 12 | 500 | 5208.819 |
| 12 | 100 | 5232.009 |
| 8 | 500 | 5259.124 |
| 11 | 100 | 5290.158 |
| 10 | 100 | 5312.424 |
Corremos un nuevo modelo utilizando los datos del mejor de los bosques aleatorios obtenido.
indice_forests <- which.min(modelos_forests$MSE_Valid)
set.seed(576327)
bosque_mejor <- randomForest(cnt ~ .,
data = train_sin,
ntree = modelos_forests[indice_forests,2],
mtry = modelos_forests[indice_forests,1],
importance=TRUE)
bosque_mejor##
## Call:
## randomForest(formula = cnt ~ ., data = train_sin, ntree = modelos_forests[indice_forests, 2], mtry = modelos_forests[indice_forests, 1], importance = TRUE)
## Type of random forest: regression
## Number of trees: 500
## No. of variables tried at each split: 10
##
## Mean of squared residuals: 1467.839
## % Var explained: 95.35Habiendo obtenido el mejor modelo de bosques aleatorios, ahora podemos ver cómo se desempeña. Para lograr lo anterior, se elaboran tres gráficas que aparecen a continuación. La primera de ella parece mostrarnos que se logró un buen ajuste para los datos de valdiación. En la segunda podemos observar los residuales y podemos percatarnos de que entre el 15 y el 30 de noviembre hay un periodo de tiempo en el que el modelo subestima la cantidad de usuarios que utilizan las bicicletas. Finalmente, en la tercera gráfica podemos observar que el modelo tiende a subestimar el número de usuarios cuando el número real de estos es mayor a 500.
preds_bosque <- predict(bosque_mejor, newdata = valid_sin)
ggplot() +
geom_line(aes(x = fechas_valid$datetime, y = valid_sin$cnt), color = "blue", alpha = 0.5, size = 0.5) +
geom_line(aes(x = fechas_valid$datetime, y = preds_bosque), color = "red", alpha = 0.5, size = 0.5) +
ggtitle("Predicciones de bosques aleatorios y valores reales a través del tiempo") +
ylab("Número de usuarios") +
xlab("Fecha") +
theme_tufte(base_size = 14, base_family = "sans")
ggplot() +
geom_point(aes(x = fechas_valid$datetime, y = valid$cnt - preds_bosque)) +
geom_abline(intercept = 0, slope = 0, color = "red", size = 1.2) +
ylab("Residuales") +
xlab("Fecha") +
ggtitle("Residuales de bosques aleatorios") +
theme_tufte(base_size = 14, base_family = "sans")
ggplot() +
geom_point(aes(x = valid$cnt, y = preds_bosque)) +
geom_abline(intercept = 0, slope = 1, color = "red", size = 1.2) +
labs(x = "Usuarios", y = "Usuarios predecidos") +
ggtitle("Predicciones de bosques aleatorios vs. Valores reales") +
theme_tufte(base_size = 14, base_family = "sans")
Gradient Boosting Machine
Se modela el problema con GBM probando distintos valores para los hiperparámetros n.trees, interaction.depth, shrinkage y bag.fraction. En particular se consideran los siguientes valores de los hiperparámetros:
- n.trees: 100, 250, 500, 750, 1000
- interaction.depth: 1, 2, 3, 4
- shrinkage: 0.001,0.005, 0.01, 0.05, 0.1, 0.5, 1
- bag.fraction: 0.25, 0.5, 0.75, 1
Dado que gbm utiliza las primeras train.fraction observaciones para entrenamiento y las siguientes 1-train.fraction observaciones para hacer la validación, debemos ingresar los datos como una matriz de la forma \(data=[train, valid]'\) y especificar el porcentaje que ocupa la muestra de entrenamiento.
#Función gbm
ajustar_boosting<-function(data.frame){
out_fun<-function(n.trees=500, interaction.depth=1,
shrinkage=0.1, bag.fraction=1,...){
set.seed(576327)
mod_boosting<-gbm(cnt~.,
data=data.frame,
distribution='gaussian', # Reg. logística
n.trees=n.trees, #no. arboles
interaction.depth =interaction.depth, #Grado de interacción variables
shrinkage=shrinkage, # Tasa de aprendizaje
bag.fraction=bag.fraction, # % de datos de entrenamiento
train.fraction = nrow(train_sin)/nrow(data.frame), # Procentaje para entrenamiento
keep.data=TRUE)
mod_boosting
}
out_fun
}
# Evaluación Modelos
Eval_gbm<-function(data){
fun_eval<-function(gbm_ajust){
# se seleccionan las variables indpe.
data_x<-data%>% select(-cnt)
# Se obtienen las prediciones
cnt_hat<-predict.gbm(gbm_ajust,newdata=data_x,n.trees=gbm_ajust$n.trees,type='response')
# Se obtiene el ECM
sqrt(mean((data$cnt-cnt_hat)^2))
}
}En total se corren 560 modelos (correspondientes a todas las posibles coibnaciones de hiperparámetros) y se elige el modelo que tenga menor error de validación. La siguiente tabla muestra los resultados obtenidos.
# Se obtienen todas las combinaciones posibles de hiperparámetros
# params<-list(n.trees=c(100,250,500,750,1000),
# interaction.depth=c(1,2,3,4),
# shrinkage=c(0.001,0.005,0.01,0.05,0.1,0.5,1),
# bag.fraction=c(0.25,0.5,0.75,1))%>%expand.grid
params<-list(n.trees=c(100,250,500,1000),
interaction.depth=c(1,2,3),
shrinkage=c(0.001,0.01,0.1,1),
bag.fraction=c(0.25,0.5,0.75,1))%>%expand.grid
# params<-list(n.trees=c(100,250),
# interaction.depth=c(1),
# shrinkage=c(0.1),
# bag.fraction=c(0.25))%>%expand.grid
# Se juntan las muestras de entrenamiento y validación.
data_train_gbm<-rbind(train,valid)
# %%%%%%%%%%%%%%%%%%%%%% Calibración hiperparámetros modelo cocina %%%%%%%%%%%%%%%%%%%%%%
# Se evalúan la funciones ajustar_boosting en los datos de entrenamiento y validacion
gbm_hiper<-ajustar_boosting(data_train_gbm)
# Se evalua la función Eval_gbm en la muestra de entrenamiento y de validación.
Eval_train_gbm<-Eval_gbm(train)
Eval_valid_gbm<-Eval_gbm(valid)
# Se corre un modelo gbm para cada combinación de hiperparámetros.
modelos_gbm<-params%>%
mutate(modelo=pmap(.,gbm_hiper))%>%
mutate(Err_Train=unlist(lapply(modelo,Eval_train_gbm)))%>%
mutate(Err_Valid=unlist(lapply(modelo,Eval_valid_gbm)))%>%
select(n.trees,interaction.depth,shrinkage,bag.fraction,Err_Train,Err_Valid)%>%
arrange(Err_Valid)
# Se imprime la tabla de modelos.
modelos_gbm %>% head(10) %>% knitr::kable(digits = 2)| n.trees | interaction.depth | shrinkage | bag.fraction | Err_Train | Err_Valid |
|---|---|---|---|---|---|
| 1000 | 3 | 1 | 1.00 | 31.44 | 68.66 |
| 500 | 3 | 1 | 1.00 | 37.06 | 70.36 |
| 500 | 3 | 1 | 0.75 | 38.38 | 71.53 |
| 250 | 3 | 1 | 0.75 | 43.93 | 72.15 |
| 1000 | 2 | 1 | 1.00 | 43.48 | 72.39 |
| 1000 | 3 | 1 | 0.75 | 34.00 | 72.45 |
| 1000 | 2 | 1 | 0.75 | 44.77 | 74.26 |
| 500 | 2 | 1 | 1.00 | 46.91 | 74.28 |
| 250 | 3 | 1 | 1.00 | 43.87 | 74.91 |
| 1000 | 2 | 1 | 0.50 | 48.31 | 75.05 |
Se toman los hiperparámetros de el “mejor modelo” bajo el método gbm y se obtienen los siguientes resultados.
# Se corre el modelo gbm con los hiperparámetros que dan el mejor modelo
mod_gbm_best<-gbm_hiper(n.trees=modelos_gbm$n.trees[1],
interaction.depth=modelos_gbm$interaction.depth[1],
shrinkage=modelos_gbm$shrinkage[1],
bag.fraction=modelos_gbm$bag.fraction[1])
# Se obtiene el desempeño del modelo en cada árbol agregado.
graf_eval<-data.frame(train=sqrt(mod_gbm_best$train.error),
valid=sqrt(mod_gbm_best$valid.error),
n_arbol=1:length(mod_gbm_best$train.error))%>%
gather(tipo,error,-n_arbol)
# Se grafica el desemeño del modelo.
ggplot(graf_eval,aes(x=n_arbol,y=error,colour=tipo,group=tipo))+
geom_line() +
theme_tufte()
Se obtiene el ajuste del mejor modelo.
# Se obtienen las predicciones de entrenamiento y prueba
cnt_Train_gbm<-predict.gbm(mod_gbm_best,newdata=train%>%select(-cnt),n.trees=mod_gbm_best$n.trees,type='response')
cnt_Valid_gbm<-predict.gbm(mod_gbm_best,newdata=valid%>%select(-cnt),n.trees=mod_gbm_best$n.trees,type='response')
# Se grafican los datos .
data_ajuste_gbm<-data.frame(Muestra=c(rep("Train",nrow(train)),rep("Valid",nrow(valid))),
Fecha=fechas_train %>% rbind(fechas_valid),
cnt_Obs=c(train$cnt,valid$cnt),
cnt_Ajust=c(cnt_Train_gbm,cnt_Valid_gbm))
ggplot() +
geom_line(aes(x = fechas_valid$datetime, y = valid_sin$cnt), color = "blue", alpha = 0.5, size = 0.5) +
geom_line(aes(x = fechas_valid$datetime, y = cnt_Valid_gbm), color = "red", alpha = 0.5, size = 0.5) +
ggtitle("Predicciones de bosques aleatorios y valores reales a través del tiempo") +
theme_tufte(base_size = 14, base_family = "sans")
ggplot() +
geom_point(aes(x = fechas_valid$datetime, y = valid_sin$cnt - cnt_Valid_gbm)) +
geom_abline(intercept = 0, slope = 0, color = "red", size = 1.2) +
ggtitle("Residuales de bosques aleatorios") +
theme_tufte(base_size = 14, base_family = "sans")
ggplot(data_ajuste_gbm,aes(x=cnt_Obs,y=cnt_Ajust,colour=Muestra)) +
geom_point() +
geom_abline(slope=1,intercept=0,lwd=1.5) +
theme_tufte()
Comparación de modelos
En esta sección se lleva a cabo la comparación final de los modelos utilizando la muestra de prueba. El modelo que posea el menor error cuadrático medio de prueba será el que se seleccionará como el mejor modelo.
Regresión Lineal
Calculamos el error cuadrático medio del modelo lineal.
preds_prueba_lineal <- predict(modelo_lineal, newdata = test_sin)
resids_prueba_lineal <- test_sin$cnt - preds_prueba_lineal
MSEp_lineal <- (resids_prueba_lineal**2) %>% mean()
MSEp_lineal## [1] 11706.01Ahora podemos ver cómo ajusta el modelo a los datos de prueba.
ggplot() +
ggtitle("Ajuste del modelo lineal") +
ylab("Número de usuarios") +
xlab("Fecha") +
geom_line(aes(x = fechas_test$datetime, y = modelos_lass$y_test, colour = "Observado")) +
geom_line(aes(x = fechas_test$datetime, y = preds_prueba_lineal, colour = "Ajustado")) +
scale_colour_manual("",
breaks = c("Observado", "Ajustado"),
values = c("blue", "red")) +
theme_tufte(base_size = 14, base_family = "sans")
Regresión Lineal con regularización
Calculamos el error cuadrático medio del modelo lineal con regularización. En este caso la regularización afectó de manera negativa al modelo.
preds_prueba_lasso <- predict(modelos_lass$mejor, newx = modelos_lass$x_test)
resids_prueba_lasso <- modelos_lass$y_test - preds_prueba_lasso
MSEp_lasso <- (resids_prueba_lasso**2) %>% mean()
MSEp_lasso## [1] 11708.74Ahora vemos cómo ajusta el modelo a los datos de prueba.
ggplot() +
ggtitle("Ajuste del modelo LASSO") +
ylab("Número de usuarios") +
xlab("Fecha") +
geom_line(aes(x = fechas_test$datetime, y = modelos_lass$y_test, colour = "Observado")) +
geom_line(aes(x = fechas_test$datetime, y = preds_prueba_lasso, colour = "Ajustado")) +
scale_colour_manual("",
breaks = c("Observado", "Ajustado"),
values = c("blue", "red")) +
theme_tufte(base_size = 14, base_family = "sans")
Red Neuronal
Calculamos el error cuadrático medio del modelo de redes neuronales.
preds_prueba_neuronales <- modelos_neuronales$models[[9]] %>% predict_on_batch(modelos_lass$x_test)
resids_prueba_neuronales <- modelos_lass$y_test - preds_prueba_neuronales
MSEp_neuronales <- (resids_prueba_neuronales**2) %>% mean()
MSEp_neuronales## [1] 4569.598Ahora observamos cómo ajusta el modelo a los datos de prueba.
ggplot() +
ggtitle("Ajuste de la red neuronal") +
ylab("Número de usuarios") +
xlab("Fecha") +
geom_line(aes(x = fechas_test$datetime, y = modelos_lass$y_test, colour = "Observado")) +
geom_line(aes(x = fechas_test$datetime, y = preds_prueba_neuronales, colour = "Ajustado")) +
scale_colour_manual("",
breaks = c("Observado", "Ajustado"),
values = c("blue", "red")) +
theme_tufte(base_size = 14, base_family = "sans")
Bosques Aleatorios
Calculamos el error cuadrático medio del modelo de bosques aleatorios.
preds_prueba_bosques <- predict(bosque_mejor, newdata = test_sin)
resids_prueba_bosques <- test_sin$cnt - preds_prueba_bosques
MSEp_bosques <- resids_prueba_bosques**2 %>% mean()
MSEp_bosques## [1] 5190.826Ahora podemos ver cómo ajusta el modelo a los datos de prueba.
ggplot() +
ggtitle("Ajuste de bosques aleatorios") +
ylab("Número de usuarios") +
xlab("Fecha") +
geom_line(aes(x = fechas_test$datetime, y = modelos_lass$y_test, colour = "Observado")) +
geom_line(aes(x = fechas_test$datetime, y = preds_prueba_bosques, colour = "Ajustado")) +
scale_colour_manual("",
breaks = c("Observado", "Ajustado"),
values = c("blue", "red")) +
theme_tufte(base_size = 14, base_family = "sans")
Gradient Boosting Machine
Se obtiene el ajuste del mejor modelo gbm para los datos de prueba.
# Se evalua la función Eval_gbm en la muestra de prueba
Eval_test_gbm<-Eval_gbm(test)
# Se obtienen las predicciones de cnt para la muestra de prueba
cnt_Test_gbm<-predict.gbm(mod_gbm_best,newdata=test %>% select(-cnt),n.trees=mod_gbm_best$n.trees, type='response')
# Se obtiene el error de prueba
TestErr_gbm<-Eval_test_gbm(mod_gbm_best)
TestErr_gbm**2## [1] 4658.615Ahora observamos cómo ajusta el modelo a los datos de prueba.
# Se grafican las observaciones
data_Pred_gbm<-data.frame(Fecha=fechas_test,
cnt_Obs=test$cnt,
cnt_Ajust=cnt_Test_gbm)
ggplot(data_Pred_gbm,aes(x=cnt_Obs,y=cnt_Ajust))+geom_point(col='olivedrab')+
geom_abline(slope=1,intercept=0,lwd=1.5) +
theme_tufte(base_size = 14, base_family = "sans")
ggplot() +
ggtitle("Ajuste de Gradient Boosting Machine") +
ylab("Número de usuarios") +
xlab("Fecha") +
geom_line(aes(x = fechas_test$datetime, y = modelos_lass$y_test, colour = "Observado")) +
geom_line(aes(x = fechas_test$datetime, y = cnt_Test_gbm, colour = "Ajustado")) +
scale_colour_manual("",
breaks = c("Observado", "Ajustado"),
values = c("blue", "red")) +
theme_tufte(base_size = 14, base_family = "sans")
Conclusión
En la siguiente se pueden observar los errores cuadráticos medios de prueba de los modelos anteriores. La red neuronal es el modelo con el menor error cuadrático medio de prueba, por lo tanto es el modelo que ha sido seleccionado por el equipo. Además, en la sección anterior se puede ver que las predicciones de este modelo son las que mejor se ajustan a los datos observados. A pesar de esto, el modelo falla al predecir el número de usuarios en los últimos días de diciembre. Esto se puede deber a que en ese periodo la gente suele utilizar menos la bicicleta por encontrarse en vacaciones navideñas. En estudios posteriores que cuenten con información de más años se podría incorporar esta variable, independientemende de la variable holiday, para verificar que esa es la razón por la que el modelo no predice correctamente los valores de esas fechas.
data.frame(Modelo = c("Lineal", "LASSO", "Red Neuronal", "Bosque Aleatorio", "GBM"),
`MSE prueba` = c(MSEp_lineal, MSEp_lasso, MSEp_neuronales,
MSEp_bosques, TestErr_gbm**2)) %>%
knitr::kable(digits = 0)| Modelo | MSE.prueba |
|---|---|
| Lineal | 11706 |
| LASSO | 11709 |
| Red Neuronal | 4570 |
| Bosque Aleatorio | 5191 |
| GBM | 4659 |